Skip to content

Conversation

@fsicignano
Copy link
Contributor

New Hubbard U term implementation.

Key features: - New hubbard.jl program
- New operator added in operators.jl
- OrbitalManifold struct has been moved to hubbard.jl
- self_consistent_field.jl has been updated to store the Hubbard n matrix and the occupation
- Hubbard term has been added to the tests in hamiltonian_consistency
- The pseudopotential used in hamiltonian_consistency has been changed with one having atomic orbitals

Copy link
Collaborator

@Technici4n Technici4n left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't looked at the math yet, but have some comments already

test_consistency_term(ExternalFromFourier(X -> abs(norm(X))))
test_consistency_term(LocalNonlinearity-> ρ^2))
test_consistency_term(Hartree())
test_consistency_term(Hubbard(DFTK.OrbitalManifold(;species=:Si, label="3S"), 1.0))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
test_consistency_term(Hubbard(DFTK.OrbitalManifold(;species=:Si, label="3S"), 1.0))
test_consistency_term(Hubbard(DFTK.OrbitalManifold(; label="3D"), 1.0))

Spacing, and let's use P to test multiple m quantum numbers. And does it work if species is removed as well?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No it doesn't work without species, should it? The intended way should be (in the multiU update that we'll make after this pr is resolved) to add each of the manifolds one by one (so to get a 3d manifold you explicitly ask for all of the single atomic manifolds having 3d orbitals)

- New handling of measure unit for U value;
- Hubbard matrix has been renamed everywhere as "nhubbard";
- Function "compute_hubbard_nIJ" has been renamed as "compute_nhubbard";
- New documentation for "compute_nhubbard" function;
- New documentation for "Hubbard" struct;
- New test on Wigner matrices added in tests/hubbard.jl;
- Refactoring inside apply! for HubbardUOperator;
- Several lines have been shortened for compactness;
- Unutilized functions have been removed from terms/hubbard.jl;
- Function reshape_hubbard_proj has been updated.
    - New handling of measure unit for U value;
    - Hubbard matrix has been renamed everywhere as "nhubbard";
    - Function "compute_hubbard_nIJ" has been renamed as "compute_nhubbard";
    - New documentation for "compute_nhubbard" function;
    - New documentation for "Hubbard" struct;
    - New test on Wigner matrices added in tests/hubbard.jl;
    - Refactoring inside apply! for HubbardUOperator;
    - Several lines have been shortened for compactness;
    - Unutilized functions have been removed from terms/hubbard.jl;
    - Function reshape_hubbard_proj has been updated.
U :: Float64
end
function (hubbard::Hubbard)(basis::AbstractBasis)
isempty(hubbard.U) && return TermNoop()
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's weird. And it's probably not what you're looking for:

julia> isempty(0)
false

@fsicignano fsicignano force-pushed the hubbard-v2 branch 2 times, most recently from 15ea5c0 to 053450f Compare September 30, 2025 14:31
- Kwarg use_nlcc=false has been added for all meta_gga tests;
- Kwarg spin_polarization=:none added to anyonic test, to avoid failure
  from @Assert length(basis.kpoints) == 1;
- Tolerance for condition number in wigner_d_matrix increased to 100,
  neq decreased to 2*l+2.
@mfherbst mfherbst changed the title Hubbard v2 Add Hubbard corrections and DFT+U Oct 7, 2025
natoms = length(manifold_atoms) # Number of atoms of the species in the manifold
nhubbard = Array{Matrix{Complex{T}}}(undef, nspins, natoms, natoms)
p_I = [Vector{Matrix{Complex{T}}}(undef, natoms) for ik in 1:length(basis.kpoints)]
# Very low-level, but works
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Extract l from first label, make sure that all the labels have the same l, and that each block of 2l+1 consecutive labels is the same except for m.

@Technici4n
Copy link
Collaborator

@antoine-levitt you said you wanted to have a look, PR is ready!

@Technici4n
Copy link
Collaborator

@antoine-levitt do you still want to have a look or should we go ahead with this PR?

@antoine-levitt
Copy link
Member

Sorry I didn't see your previous message (got lost in the flood of notifications). I'll take a look!

Copy link
Member

@antoine-levitt antoine-levitt left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks nice! I just have an issue with the manifold selection API. Needs to be clearer that one atomic orbital = one Hubbard term = one l. I would suggest just

struct Hubbard
U
element::Element
atomic_orbital::String
end
struct TermHubbard
U
P
end

and getting rid of the manifold thing entirely. If you want to +U only selected atoms, then just make them different element types.

stating whether the orbital belongs to the manifold.
"""
@kwdef struct OrbitalManifold
iatom ::Union{Int64, Nothing} = nothing
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this just an int? If I want to put a +U on two atom types, should I have two hubbard terms?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At the moment it is not possible to have more than one +U correction, but I have a different version of the code which allows "multi-U" calculations with different atom types. I was planning to add it as a different PR once this one is closed.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think since it very much determines the structure of the code we should do it now

where n or ibnd is the band index, ``weights[ik ibnd] = kweights[ik] * occupation[ik, ibnd]``
and ``Pᵢₘ₁`` is the pseudoatomic orbital projector for atom i and orbital m₁
(usually just the magnetic quantum number, since l is usually fixed).
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not a "usually" it's enforced in the code, right? So it's actually enforced that each Hubbard term only deals with one type of atomic orbital. This needs to be enforced by the design.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the code this is enforced, yes, because the symmetrization of the nhubbard is l-dipendent.

Copy link
Member

@mfherbst mfherbst left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @fsicignano for updating this so quickly. I think it indeed makes for a cleaner code.

@Technici4n
Copy link
Collaborator

There is one more thing: currently the example takes a few minutes to run, maybe its discretization should be made coarser.

@mfherbst
Copy link
Member

There is one more thing: currently the example takes a few minutes to run, maybe its discretization should be made coarser.

Yes, I agree. If we can, we should. But with magnetism that's sometimes not possible.

@Technici4n
Copy link
Collaborator

I think we are good now!

@Technici4n
Copy link
Collaborator

@antoine-levitt now is the time if you want to have another look :)

@antoine-levitt
Copy link
Member

LGTM, thanks for doing this!

@mfherbst mfherbst merged commit 908e3c8 into JuliaMolSim:master Nov 14, 2025
8 of 10 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants